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Abstract 

In  this  work,  a  three-dimensional,  steady-state,  multi- agglomerate  model  of  cathode  catalyst  layer  in  polymer  electrolyte  membrane  (PEM)  fuel 
cells  has  been  developed  to  assess  the  activation  polarization  and  the  current  densities  in  the  cathode  catalyst  layer.  A  finite  element  technique  is 
used  for  the  numerical  solution  to  the  model  developed.  The  cathode  activation  overpotentials,  and  the  membrane  and  solid  phase  current  densities 
are  calculated  for  different  operating  conditions.  Three  different  configurations  of  agglomerate  arrangements  are  considered,  an  in-line  and  two 
staggered  arrangements.  All  the  three  arrangements  are  simulated  for  typical  operating  conditions  inside  the  PEM  fuel  cell  in  order  to  investigate 
the  oxygen  transport  process  through  the  cathode  catalyst  layer,  and  its  impact  on  the  activation  polarization.  A  comprehensive  validation  with  the 
well-established  two-dimensional  “axi- symmetric  model”  has  been  performed  to  validate  the  three-dimensional  numerical  model  results.  Present 
results  show  a  lowest  activation  overpotential  when  the  agglomerate  arrangement  is  in-line.  For  more  realistic  scenarios,  staggered  arrangements,  the 
activation  overpotentials  are  higher  due  to  the  slower  oxygen  transport  and  lesser  passage  or  void  region  available  around  the  individual  agglomerate. 
The  present  study  elucidates  that  the  cathode  overpotential  reduction  is  possible  through  the  changing  of  agglomerate  arrangements.  Hence,  the 
approaches  to  cathode  overpotential  reduction  through  the  optimization  of  agglomerate  arrangement  will  be  helpful  for  the  next  generation  fuel 
cell  design. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Polymer  electrolyte  membrane  (PEM)  fuel  cell  is  a  promis¬ 
ing  candidate  for  the  next  generation  power  sources  due  to  its 
high  power  density,  low  operating  temperature,  quick  start-up, 
and  fast  dynamic  response.  Most  importantly,  its  zero  emis¬ 
sion  capabilities  open  up  opportunities  to  wide  practical  use  in 
portable,  mobile,  and  stationary  cogeneration  applications  [1]. 
Even  though  substantial  improvements  have  been  made  over 
the  past  few  years  in  the  cell  design  and  material  utilization  in 
PEM  fuel  cells,  several  technical  barriers  still  exist  that  prevent 
the  PEM  fuel  cells  from  commercialization.  Among  the  various 
obstacles,  the  most  important  are  low  cell  performance  due  to 
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high  polarization,  and  high  cost  due  to  platinum  (Pt)  catalyst 
used  in  the  catalyst  layer  (CL).  Hence,  there  has  been  consider¬ 
able  interest  in  the  modeling  and  optimization  of  PEM  fuel  cells 
aiming  at  performance  improvement  and  cost  reduction  [2-17]. 

The  performance  of  a  PEM  fuel  cell  is  mainly  dictated 
by  ohmic,  activation,  and  concentration  overpotentials.  At  the 
most  common  operating  ranges,  ohmic  and  activation  over¬ 
potentials  are  dominant  over  the  concentration  overpotential, 
or  simply  existence  of  concentration  overpotential  is  negli¬ 
gible.  The  estimation  of  the  ohmic  overpotential  has  already 
been  well-established  that  can  be  determined  from  the  avail¬ 
able  experimental  data  or  from  the  empirical  relation  for  the  cell 
polarization  curve  of  PEM  fuel  cells  [6-8,11,18].  On  the  other 
hand,  the  activation  overpotential  shows  more  complex  nature 
due  to  its  dependency  on  the  catalyst  layer  structure;  for  instance, 
whether  it  has  agglomerate  or  macro-homogeneous  structure, 
composition  of  the  catalyst  layer,  types  of  the  catalysts  used, 
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and  how  the  reactants  transport  in  the  catalyst  layer.  Therefore, 
activation  overpotential  cannot  be  estimated  easily.  Further¬ 
more,  activation  overpotential  in  the  anode  catalyst  layer  is  very 
small  compared  to  the  activation  overpotential  in  the  cathode 
catalyst  layer;  hence,  the  activation  overpotential  in  the  cath¬ 
ode  catalyst  layer  has  major  influence  on  the  cell  performance. 
Among  the  various  catalyst  layer  models  investigated,  three  dif¬ 
ferent  models  of  cathode  catalyst  layer  have  been  established 
in  the  past  decade,  namely,  thin-film,  macro-homogeneous,  and 
agglomerate  model.  In  the  thin-film  model,  the  catalyst  particles 
are  embedded  on  the  thin-film  of  polymer  membrane  [19,20]; 
whereas  in  the  macro-homogeneous  model,  the  cathode  catalyst 
layer  is  considered  as  a  homogeneous  matrix  of  supported  cat¬ 
alyst  platinum,  polymer  electrolyte,  and  void  spaces  [4,5,9,10]. 
In  the  agglomerate  model,  the  catalyst  layer  is  considered  as  a 
uniform  matrix  of  catalyst  agglomerates,  which  is  surrounded  by 
the  gas  pores.  Each  of  these  catalyst  agglomerates  are  assumed 
to  be  homogeneous  mixture  of  catalysts,  polymer  electrolytes, 
and  void  spaces  as  well  [15,21-27].  In  addition,  there  are  several 
other  models  that  have  also  been  developed,  namely,  cylindrical 
agglomerate  [28],  ordered  catalyst  layer  [29],  and  non-uniform 
catalyst  layer  [30].  It  should  be  noted  here  that  experimental 
studies  showed  that  the  agglomerate  model  might  be  a  close 
approximation  to  model  the  catalyst  layer  for  PEM  fuel  cells 
[22,31]. 

On  the  other  hand,  in  terms  of  the  computational  efficiency  in 
numerical  modeling,  three  different  approaches  to  catalyst  layer 
modeling  would  be  more  viable.  In  a  most  simplified  approach, 
the  catalyst  layer  is  considered  as  an  ultra-thin  layer  between 
the  membrane  and  gas  diffusion  layer  (GDL),  ignoring  the  reac¬ 
tion  kinetics  and  transport  processes  within  the  catalyst  layer 
[32].  This  is  useful,  particularly  for  the  three-dimensional  full 
cell  model  due  to  the  limited  computing  power.  Hence,  the  cat¬ 
alyst  layer  is  assumed  as  a  source  or  sink  boundary  condition 
that  is  simple  to  implement  for  predicting  the  typical  polarization 
curve,  optimizing  the  cell  design  parameters,  and  operating  con¬ 
ditions.  In  the  second  approach,  the  catalyst  layer  is  considered 
as  more  realistic  thin  layer  of  catalyst  particles  and  electrolyte 
membrane  or  the  agglomerate  of  catalyst  particles  and  elec¬ 
trolyte  membrane  sandwiched  between  the  membrane  and  gas 
diffusion  layer  [4,5,9,10,15,23-25,30,33].  All  of  these  studies 
used  the  one-dimensional  approach  for  the  catalyst  layer  or  for 
the  individual  catalyst  agglomerate.  Recently,  two-dimensional 
model  results  are  presented  for  the  agglomerate  model  [27]. 
Once  again,  the  individual  agglomerate  is  accounted  by  consid¬ 
ering  diffusion  in  the  radial  direction  only.  However,  to  improve 
the  performance  of  the  cathode  catalyst  layer  and  reduce  the  cost 
associated  with  the  Pt-catalyst,  it  is  required  to  study  transport 
processes  in  the  cathode  catalyst  layer.  Furthermore,  transport 
processes  in  the  cathode  catalyst  layer  are  largely  depend  on 
the  structures  of  the  catalyst  layer  or  the  arrangements  of  cat¬ 
alyst  agglomerates  in  the  cathode  catalyst  layer.  Therefore,  the 
third  and  most  practical  approach  would  be  detailed  modeling 
of  the  cathode  catalyst  layer  using  three-dimensional  approach. 
None  of  the  previous  studies  has  reported  the  three-dimensional 
nature  of  the  agglomerate  catalyst  layer  that  would  be  an  accurate 
approximation  of  a  practical  catalyst  layer. 


In  this  study,  a  three-dimensional  (3D)  agglomerate  model  of 
cathode  catalyst  layer  in  a  PEM  fuel  cell  is  developed  to  study 
the  activation  overpotential  in  the  cathode  catalyst  layer.  The 
effect  of  agglomerate  arrangements  on  the  activation  overpoten¬ 
tial  of  PEM  fuel  cells  has  been  investigated  for  three  different 
types  of  agglomerate  arrangements,  namely,  in-line  agglomerate 
arrangement  as  Case-I,  and  two  staggered  agglomerate  arrange¬ 
ments  as  Case-II  and  Case-Ill.  The  catalyst  layer  geometry  is 
generated  assuming  that  the  agglomerates  are  aligned  along  the 
thickness  of  the  catalyst  layer  in  the  first  case  and  then  by  consid¬ 
ering  staggered  arrangements  in  the  subsequent  cases.  Since  the 
governing  equations  are  nonlinear  partial  differential  equations 
and  the  catalyst  layer  has  a  complex  geometry,  finite  element 
technique  is  used  to  solve  the  governing  equations  due  to  its 
ability  to  handle  complex  geometrical  domains  [34,35].  The 
simulation  results  presented  here  show  a  considerable  insight 
on  how  the  activation  overpotential  changes  with  the  arrange¬ 
ment  of  catalyst  agglomerates  that  will  eventually  be  helpful  for 
the  better  understanding  of  PEM  fuel  cell  performance  and  its 
design. 

2.  Model  description 

A  typical  PEM  fuel  cell  is  considered  that  consists  of  a 
cathode  and  an  anode  electrode  with  a  proton  conducting  mem¬ 
brane  as  the  electrolyte  sandwiched  in  between.  Generally,  the 
thickness  of  the  electrodes  and  membrane  are  approximately 
250  pan.  Each  of  these  electrodes  also  consists  of  approximately 
10  [Jim  (or  thinner)  catalyst  layer  between  the  electrode  and 
the  membrane,  known  as  the  anode  catalyst  layer  and  cathode 
catalyst  layer,  respectively.  Typically,  humidified  H2  gas  is  sup¬ 
plied  under  pressure  into  the  anode  flow  channel  which  diffuses 
through  the  porous  electrode  until  it  reaches  the  anode  catalyst 
layer  and  forms  protons  (H+)  and  electrons  via  electro-oxidation 
reaction  at  the  catalyst  surface.  The  protons  are  transferred 
through  the  membrane  to  the  cathode  catalyst  layer,  and  the 
electrons  are  transported  via  the  external  circuit  to  the  cathode. 
Conversely,  humidified  O2  gas  or  air  is  supplied  to  the  cathode 
flow  channel  where  O2  gas  diffuses  through  the  porous  electrode 
until  it  reaches  the  cathode  catalyst  layer  and  forms  water  react¬ 
ing  with  protons  and  electrons.  The  overall  electro-chemical 
reaction  occurring  in  the  PEM  fuel  cell  can  be  represented  by 
the  following  reaction: 

H2  +  ^C>2  ->  H2O  +  Heat  +  Electric  energy  (1) 

2.7.  Governing  equations 

In  this  study,  a  three-dimensional  modeling  domain  is  consid¬ 
ered,  which  is  geometrically  identical  to  our  previous  study  [18], 
except  the  catalyst  layer  structure.  A  schematic  diagram  of  the 
cathode  catalyst  layer  under  consideration  is  shown  in  Fig.  1. 
Here,  the  catalyst  agglomerates  are  considered  as  a  homoge¬ 
neous  mixture  of  electrolyte  membrane,  supported  Pt-catalyst, 
and  void  space.  The  volume  fractions  of  these  components  can 
be  varied  as  can  the  effective  surface  area  of  catalyst  that  can 
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be  characterized  by  different  loadings  and  catalyst  types.  The 
overall  reaction  in  the  catalyst  agglomerate  is  taken  as 

02  +  4H+  +  4e“  -*  2H20  (2) 

Assuming  the  cell  is  operating  in  steady-state  condition  and 
the  membrane  is  fully  humidified,  the  conservation  equation  for 
oxygen  in  the  cathode  catalyst  layer  (inside  the  agglomerates) 
can  be  written  as 

V  •  (d^VCo2)  +  ^orr  =  o  (3) 

where  D is  the  effective  diffusion  coefficient  of  the  trans¬ 
ported  oxygen  that  is  calculated  using  Bruggeman  correction 
from  the  bulk  diffusion  coefficient,  Do2 ,  and  the  corresponding 
void  fraction,  0;  as  [36] 

ml = f,2D<h  (4) 

Since  there  is  no  reaction  outside  the  agglomerates,  the  oxygen 
reduction  reaction  rate  or  the  rate  of  electro-chemical  reaction 
(^orr)  is  equal  to  zero  outside  the  agglomerates  in  Eq.  (3);  there¬ 
fore,  diffusion  is  the  mechanism  for  oxygen  transport  in  the 
void  region  outside  the  agglomerates.  Conversely,  the  oxygen 
reduction  reaction  rate  inside  the  agglomerates  is  given  by  the 
Butler-Volmer  equation  as  [1] 


Here,  the  catalyst  reactive  surface  area  per  unit  volume  (A„)  is 
a  function  of  the  catalyst  mass  loading  per  unit  area  of  cathode 
(mpt),  the  catalyst  surface  area  per  unit  mass  of  the  catalyst  (As), 
and  the  thickness  of  the  catalyst  layer  (<$c),  which  is  calculated 
using  the  following  relation  [18]: 

.  AsmPt 

An,  =  -  (6) 

c 


The  reference  current  density  (/o,ref)  at  the  reference 
concentration  of  Co2,ref  is  calculated  using  the  experimental 
data  of  Parthasarathy  et  al.  [38]  and  the  reference  oxygen 
concentration,  Co2,ref>  is  taken  as  1.2  mol  m-3 [10,38].  The 
experimental  data  of  the  reference  exchange  current  density  in 
A  cm-2  for  oxygen  reduction  in  Nation®  were  correlated  with 
the  cell  temperature  in  Kelvins  by  [10] 

4001 

log10(/o,ref)  =  3.507  -  —  (7) 

In  Eq.  (5),  y  is  the  overall  reaction  order,  ofa  and  ac  are  the 
apparent  transfer  coefficients  for  the  anodic  and  cathodic 
reactions,  respectively,  and  r/nct  represents  the  activation 
overpotential  for  the  electro-chemical  reactions. 

It  is  assumed  that  the  electric  current  flowing  through  the 
catalyst  layer  obeys  the  Ohm’s  law,  which  relates  the  electrical 
current  density  to  the  potential  gradient  as 

Ji  =  -ofVVt  (8) 

where  crff  is  the  effective  conductivity  of  the  medium  through 
which  the  current  travels  and  index  i  refers  the  phases  in  the 
catalyst  layer.  Since  the  catalyst  agglomerate  in  the  catalyst 
layer  consists  both  the  electrolyte  membrane  phase  and  the  solid 
catalyst,  Eq.  (8)  can  be  written  for  the  catalyst  agglomerates  as 

Jm  =  —  for  membrane  phase  current  (9) 

Js  = —crffWs  for  solid  phase  current  (10) 

where  Jm  and  Js  are  the  membrane  phase  current  density  and  the 
solid  phase  current  density,  and  crff  are  the  effective  protonic 
conductivity  of  the  membrane  and  the  effective  electronic  con¬ 
ductivity  of  the  carbon  support  in  the  agglomerate,  and  &s 
are  the  membrane  potential  and  the  solid  phase  potential,  respec¬ 
tively.  In  addition,  the  current  conservation  equation  yields  [37] 

V/m  +  V/s  =  0  (11) 

In  the  cathode  catalyst  layer,  the  protonic  current  is  defined 
as  the  flow  of  positively  charged  particles  and  the  electronic 
current  is  defined  as  the  flow  of  negatively  charged  particles. 
Since  the  flow  of  electrons  depends  on  the  flow  of  oxygen  in  the 
cathode  catalyst  layer  through  Eq.  (2),  the  relationship  between 
the  solid  phase  current  and  the  oxygen  reduction  reaction  rate 
can  be  related  as 

V  •  Js  =  -/i39t orr  (12) 

where  n  is  the  number  of  electrons  transferred  in  the  cathodic 
reaction  in  Eq.  (2). 

By  combining  Eqs.  (9)-(12),  the  membrane  phase  poten¬ 
tial  and  the  solid  phase  potential  equations  within  the  catalyst 
agglomerates  are  summarized  as  follows: 

V  •  =  —n%$i on-  for  membrane  phase  potential 

(13) 

V  •  ^<7gffV^s^  =  n^?don  for  solid  phase  potential 


(14) 
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where  and  crff  within  the  catalyst  agglomerates  are  calcu¬ 
lated  with  Bruggemann  correction  from  the  bulk  conductivity 
using  the  following  expressions: 


—  °m(/m  x  0c)  ^ 

(15) 

=  ffg(l  -  0c)3/2 

(16) 

where  lm  denotes  the  volume  fraction  of  membrane  in  the 
agglomerates,  am  and  as  are  the  bulk  conductivities  of  the  mem¬ 
brane  and  the  solid  catalyst,  respectively,  and  0C  is  the  void 
fraction  in  the  agglomerates. 

2.2.  Boundary  conditions 


n./m  =  0  for  df2  e  agglomerate-void  space  interface  (20c) 

where  Jm  is  the  membrane  current  density  at  the  membrane-CL 
interface.  Finally,  for  the  solid  phase  potential  as 

-crfs/^s  =  Js  for  d£2  e  GDL-CL  interface  (21a) 

=  0  for  df2  g  membrane-CL  interface  (21b) 

n./s  =  0  for  3 agglomerate-void  space  interface  (21c) 

where  Js  is  the  solid  phase  current  density  at  the  GDL-CL  inter¬ 
face.  Due  to  the  geometric  symmetry,  the  total  membrane  current 
density  at  the  membrane-CL  interface  is  equal  to  the  total  solid 
phase  current  density  at  the  GDL-CL  interface,  i.e., 


The  overpotential  in  the  cathode  catalyst  layer  is  defined  as 
the  potential  difference  between  the  local  value  and  the  reference 
potential.  In  order  to  be  able  to  calculate  the  overpotential  in  the 
cathode  catalyst  layer,  it  is  required  to  define  a  reference  poten¬ 
tial  in  the  catalyst  layer.  In  the  present  model,  the  solid  phase 
potential  is  considered  as  the  reference  potential  that  eventually 
simplified  the  electronic  overpotential  as  zero  and  the  protonic 
overpotential  as  the  total  activation  overpotential.  Hence,  the 
activation  potential  in  the  cathode  catalyst  layer  is  simplified  to 


*7 act  —  (12) 

Furthermore,  following  boundary  conditions  are  used  to  solve 
the  governing  partial  differential  equations.  For  the  oxygen 
transport  equation,  the  boundary  conditions  are  defined  as 

Cq2  =  Cmt  for  dfi  e  GDL-CL  interface  (18a) 

VCo2  =  0  for  3f2  e  membrane-CL  interface  (18b) 

n-(Nin  —  N0ut)  =  0 

for  3£?e  agglomerate-void  space  interface  (18c) 


where  Cmi  is  the  oxygen  concentration  at  the  GDL-CL  interface, 
3£?  represents  the  boundary  of  the  computational  domain,  and  n 
is  the  unit  normal  vector.  Nin  and  Nout  represent  the  inward 
and  outward  fluxes  at  the  agglomerate-void  space  interface, 
respectively,  that  is  defined  as 

N|  =  — AV(Co2);  (19) 

where  index  i  refers  inward  or  outward  direction.  The  calcula¬ 
tion  of  Cint  is  not  trivial  since  oxygen  concentration  is  normally 
known  in  the  flow  channel,  and  decreases  after  transporting 
through  the  GDL  to  reach  the  catalyst  layer.  However,  an  appro¬ 
priate  measure  is  required  at  the  liquid-gas  interface  for  partially 
flooded  agglomerates  in  the  catalyst  layer  since  the  reactant  gas 
is  weakly  soluble  in  liquid  water  under  the  typical  fuel  cell  oper¬ 
ating  condition,  Henry’s  law  can  be  used  to  relate  the  reactant 
concentrations  in  the  gas  and  liquid  phases.  The  details  are  given 
in  our  earlier  work  [18].  For  the  membrane  phase  potential,  the 
boundary  conditions  are 

V^m  =  0  for  3 GDL-CL  interface  (20a) 

— ^nfv^m  =  Jm  for  3f2  e  membrane-CL  interface  (20b) 


Jm  e  membrane-CL  interface  =  Js  e  GDL-CL  interface  =  Js 

(22) 

where  Js  is  the  boundary  value  of  the  current  densities  for 
the  two  interfaces  as  mentioned  above.  In  addition  to  the 
above-mentioned  boundary  conditions,  insulation  or  symmetry 
boundary  conditions  are  applied  in  the  appropriate  boundaries 
when  required. 


2.3.  Operating  conditions  and  physical  parameters 


In  the  present  investigation,  it  is  considered  that  air  is  the 
cathode  gas  and  the  concentration  of  oxygen  in  the  cathode  flow 
channel  is  uniform.  The  electrode  is  considered  dry;  hence,  the 
oxygen  diffuses  through  the  un-flooded  electrode  void  region  to 
reach  the  GDL-CL  interface.  In  addition,  the  catalyst  agglom¬ 
erates  are  considered  partially  hydrated  and  water  in  the  void 
region  around  the  agglomerates  is  considered  in  gaseous  phase 
to  simulate  the  un-flooded  scenarios  and  the  start-up  case.  The 
thickness  of  the  catalyst  layer  in  the  present  study  is  considered 
as  10  p,m  and  the  agglomerate  diameter  is  considered  as  5  [xm. 
Typically,  the  thickness  of  the  catalyst  layer  and  the  agglomerate 
size  depend  on  the  amount  of  catalyst  loading  and  the  fabrication 
methods.  It  is  found  in  our  earlier  study  that  for  the  typical  oper¬ 
ating  conditions  and  physical  parameters,  the  optimum  catalyst 
thickness  ranges  from  10  to  15  [xm  [18].  Furthermore,  scanning 
electron  micrograph  (SEM)  of  membrane  electrode  assembly 
shows  that  the  catalyst  layer  thickness  is  around  10-20  |uim  and 
the  mean  agglomerate  diameter  is  about  6  [xm  [22].  The  operat¬ 
ing  parameters  and  the  physical  properties  used  in  the  numerical 
computation  are  listed  in  Table  1 . 

It  is  worthwhile  to  note  that  the  bulk  diffusion  coefficient  of 
oxygen  is  calculated  according  to  the  following  relation  [39]: 


2^02  ,bulk  = 


_ l-*o2 _ 

(Xn2/Dq2- n2)  +  (2fH2o/A}2-H2o) 


(23) 


where  Xq2  ,  Xn2,  and  Xh2o  are  the  mole  fractions  of  oxygen, 
nitrogen,  and  water  vapor,  respectively.  The  binary  diffusion 
coefficient  of  oxygen  and  nitrogen,  Do2-n2,  is  calculated  using 
the  Chapman-Enskog  formula  and  the  binary  diffusion  coeffi¬ 
cient  of  oxygen  and  water  vapor,  Do2-h2o>  is  calculated  using 
the  Slattery-Bird  equation  [39,40].  Since  the  catalyst  agglom¬ 
erate  is  a  mixture  of  electrolyte  membrane,  solid  catalyst,  and 
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Table  1 

The  operating  and  physical  parameters  in  the  present  model  calculations  [10,18] 


Parameter  Value 


Operating  temperature,  T  (°C)  50  and  80 

Operating  pressure,  P  (atm)  1  and  3 

Electrode  thickness,  <5e  (|Jim)  250 

Catalyst  layer  thickness,  Sc  (fxm)  10 

Agglomerate  diameter  (|jim)  5 

Void  fraction  of  the  cathode  electrode,  0.4 

Fraction  of  membrane  in  the  agglomerate,  lm  0.4 

Fraction  of  flooding  in  the  agglomerate,  l\  0.5 

Catalyst  loading  per  unit  area,  rapt  (mg  cm-2)  0.2 

Fraction  of  catalyst,  /pt  0.2 

Catalyst  surface  area  per  unit  mass  of  the  catalyst,  As  (m2  g_1)  112 

Membrane  conductivity,  crm  (S  cm-1)  0.17 

Solid  catalyst  conductivity,  crs  (S  cm-1)  727 

Density  of  platinum,  ppt  (g  cm-3)  21.5 

Density  of  carbon  black,  pc  (g  cm-3)  2.0 

Anodic  transfer  coefficient,  aa  0.5 

Cathodic  transfer  coefficient,  ac  0.5 


void  region  filled  with  liquid  water  and  water  vapor,  the  overall 
diffusion  coefficient  in  the  catalyst  layer  can  be  estimated  using 
the  following  relation: 


Ad2-c  Do2-m  Ad2-h2o(i)  £)o2-H2o(g) 

where  lm,l\ ,  and  lg  are  the  fractions  of  membrane,  liquid  water, 
and  water  vapor  in  the  catalyst  agglomerate,  respectively. 

3.  Numerical  technique 

The  governing  equations,  Eqs.  (3),  (13),  and  (14)  for  inside 
and  Eq.  (3)  for  outside  of  the  agglomerates,  were  solved  sub¬ 
jected  to  the  boundary  conditions  mentioned  in  the  previous 
section.  It  should  be  pointed  out  that  two  approaches  exist  in 
literature  for  the  numerical  modeling  of  a  fuel  cell.  In  the  first 
approach,  current  densities  can  be  solved  numerically  by  speci¬ 
fying  overpotentials,  and  in  the  second  approach,  overpotentials 
are  solved  by  specifying  current  densities.  In  the  present  work, 
the  second  approach  was  used. 

The  finite  element  method  using  COMSOL  Multiphysics  ™ 
running  on  a  64-bit  Linux  CPU  with  3  GB  of  RAM  was  used. 
The  computational  domain  was  initially  discretized  into  a  tetra¬ 
hedral  mesh  and  Lagrangian  elements  of  second  order  (quadratic 


elements)  were  used.  Stationary  nonlinear  solver  was  chosen 
to  solve  the  governing  partial  differential  equations  (PDEs). 
Since  the  governing  PDEs  are  highly  nonlinear,  the  general 
form  of  solution  was  chosen  with  the  GMRES  iterative  solver 
and  the  geometric  multigrid  or  SSOR  techniques  were  used 
as  pre-conditioners.  Furthermore,  governing  equations  are  cou¬ 
pled  together  with  the  reaction  rate  term,  hence,  in  the  solution 
methodology  an  initial  solution  first  obtained  in  a  lower  mesh 
and  then  the  governing  equations  were  solved  individually  at 
the  higher  level  of  mesh.  The  solutions  were  considered  as  con¬ 
verged  solution  when  the  preset  tolerance  value  goes  below 
10-6  for  each  case.  Detail  description  of  the  solver,  the  pre¬ 
conditioners,  and  the  error  estimation  used  in  this  study  can  be 
found  in  the  COMSOL  Multiphysics  ™  user’s  guide  [41]. 

3.1.  Numerical  validation 

As  mentioned  earlier,  none  of  the  previous  studies  has  consid¬ 
ered  the  three-dimensional  agglomerate  arrangements;  hence, 
direct  comparison  is  not  possible.  Rather  a  limiting  case  of 
agglomerate  model  has  been  invoked  for  the  validation  of  the 
accuracy  of  three-dimensional  numerical  calculation,  where 
agglomerates  are  considered  in  a  cylindrical  computational 
domain.  The  advantage  of  using  such  three-dimensional  domain 
is  it  can  be  solved  as  two-dimensional  axi- symmetric  problem, 
hence,  the  accuracy  of  the  three-dimensional  calculation  can  eas¬ 
ily  be  verified  with  two-dimensional  calculation.  The  thickness 
of  the  catalyst  layer  is  chosen  as  10|mm  and  the  agglomer¬ 
ate  diameter  as  5  p,m,  i.e.,  two  agglomerates  can  exist  along 
the  thickness  of  the  catalyst  layer.  A  schematic  diagram  of 
the  three-dimensional  computational  domain  and  correspond¬ 
ing  axi-symmetric  computational  domain  used  for  the  validation 
is  shown  in  Fig.  2  along  with  the  coordinate  systems.  Here, 
the  number  of  agglomerates  was  kept  as  two,  however,  to 
maintain  sufficient  contact  between  the  agglomerates,  between 
the  agglomerates  and  GDL,  and  between  the  agglomerates 
and  membrane,  the  size  of  the  each  agglomerate  has  been 
increased  by  2%  keeping  the  centers  of  the  agglomerates  fixed. 
Boundary  conditions  for  both  the  axi-symmetric  model  and 
the  three-dimensional  model  were  identical  as  described  ear¬ 
lier. 

Fig.  3  shows  the  oxygen  concentration  profile  across  the  cata¬ 
lyst  layer  thickness  along  the  centerline  of  the  agglomerates  (line 
OO  in  Fig.  2)  for  two  current  densities  as  indicated  in  the  legend. 


GDL-CL 

Interface 


Membrane-CL 

Interface 


Fig.  2.  Schematic  of  the  computational  domain:  (a)  three-dimensional  domain,  (b)  two-dimensional  axi-symmetric  domain. 
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Fig.  3.  Oxygen  concentration  profile  in  the  cathode  catalyst  layer  along  the  x- 
axis  shown  in  Fig.  2  in  a  PEM  fuel  cell  operating  at  80  °C  and  3  atm.  Lines 
represent  the  axi-symmetric  model  results,  whereas  the  symbols  represent  the 
three-dimensional  model  results  of  the  present  investigation  for  two  current 
densities  as  indicated  in  the  legend. 

The  operating  pressure  and  temperature  of  the  fuel  cell  is  consid¬ 
ered  as  3  atm  and  80  °C,  respectively.  The  numerical  procedure 
for  the  three-dimensional  models  is  similar  as  described  earlier. 
However,  for  axi-symmetric  model,  the  advantage  of  adaptive 
mesh  refinement  technique  has  been  employed  for  better  accu¬ 
racy.  Conversely,  Fig.  4  depicts  the  activation  polarization  for 
the  three-dimensional  model  and  the  axi-symmetric  model.  As 
observed  in  Fig.  3  for  both  current  densities,  numerical  solu¬ 
tion  of  the  three-dimensional  model  shows  good  agreement 
with  the  axi-symmetric  model  results.  Furthermore,  the  acti¬ 
vation  polarization  results  also  show  excellent  agreement.  The 
accuracy  of  the  two-dimensional  finite  element  model  using 
adaptive  mesh  refinement  with  commercial  software,  like  COM- 
SOL  Multiphysics  ™,  is  well-established  [34,42] .  Furthermore, 
two-dimensional  numerical  model  required  less  number  of  grid 
to  represent  the  curve  surfaces.  Whereas,  three-dimensional 


Js  [A/cm2] 


Fig.  4.  Activation  polarization  of  a  PEM  fuel  cell  operating  at  80  °C  and  3  atm. 
Lines  represent  the  axi-symmetric  model  results,  whereas  the  symbols  represent 
the  three-dimensional  model  results  of  the  present  investigation. 


models  required  significantly  large  number  of  gird  for  proper 
representation  of  spherical  surfaces  that  is  limited  by  com¬ 
puter  memory.  This  limitation  of  computer  memory  eventually 
lower  the  accuracy  of  the  three-dimensional  numerical  cal¬ 
culation.  Since  the  comparisons  shown  here  provide  a  good 
agreement  with  each  other,  it  can  easily  be  concluded  that 
the  three-dimensional  model  results  presented  in  this  study  are 
sufficiently  accurate  for  studying  the  effect  of  catalyst  layer 
structures  on  the  performance  of  PEM  fuel  cells.  Once  again, 
these  results  show  the  accuracy  level  of  three-dimensional 
numerical  computation,  whereas  the  accuracy  of  the  mathe¬ 
matical  formulation  has  already  been  established  elsewhere 
[18]. 

3.2.  Grid  sensitivity 

In  addition  to  the  validation  shown  in  the  previous  section, 
a  comprehensive  grid  sensitivity  test  has  been  performed  for 
all  the  three  cases  (Case-I,  Case-II,  and  Case-Ill)  to  ensure  the 
results  are  independent  of  grid  sizes.  This  test  also  ensures 
the  number  of  grid  used  in  the  three-dimensional  numerical 
calculation  is  sufficient.  The  schematic  of  the  agglomerate 
arrangement  in  the  cathode  catalyst  layer  and  the  computational 
domain  for  the  three  cases  used  in  the  present  investigation  is 
shown  in  Fig.  5.  Here,  Case-I  represents  in-line  agglomerate 
arrangement,  whereas  Case-II  and  Case-Ill  depict  two  stag¬ 
gered  arrangements  (uni-directional  and  bi-directional  staggered 
arrangements).  The  orientation  of  the  catalyst  agglomerates  in 
the  catalyst  layer  for  different  agglomerate  arrangements  are 
listed  in  Table  2.  Different  level  of  meshes  have  been  con¬ 
sidered  starting  from  approximately  0.4  million  to  1  million 
degrees  of  freedom  (DOF)  for  the  three-dimensional  models. 
These  types  of  geometry  cannot  be  solved  as  axi-symmetric, 
therefore,  full  three-dimensional  computational  domain  is  used 
to  estimate  oxygen  concentrations  and  activation  overpoten¬ 
tials.  For  illustration  purpose,  the  grid  sensitivity  results  are 
shown  for  Case-I  only.  Fig.  6  illustrates  the  variation  of  the 
oxygen  concentration  and  the  activation  overpotential  in  the 
cathode  catalyst  layer  along  the  centerline  of  the  two  refer¬ 
ence  agglomerates  for  Case-I.  As  mentioned  earlier  for  Case-I, 
agglomerate  arrangements  are  in-line  in  all  directions,  hence, 
the  geometry  is  almost  identical  to  three-dimensional  domain 
shown  in  Fig.  2a,  except  the  outer  domain  is  rectangular  instead 
of  cylindrical.  All  the  simulation  parameters  are  identical  to 
those  for  Fig.  4.  Accuracy  of  the  three-dimensional  model 
has  already  been  validated  in  Figs.  3  and  4,  therefore,  these 
results  can  be  as  accurate  as  Figs.  3  and  4.  From  Fig.  6,  it 


Table  2 

Agglomerates  orientation  in  different  directions  for  the  cases  considered  in  the 
present  investigation 


x-Direction 

^-Direction 

z-Direction 

Case-I 

In-line 

In-line 

In-line 

Case-II 

In-line 

Staggered 

In-line 

Case-Ill 

In-line 

Staggered 

Staggered 
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smbrane-CL 

Interface 


Reference 

Agglomerate 


Case-I  (b) 


GDL-CL 

Interface 


Reference 

Agglomerate 


Case-I  (a) 


Reference 
Agglomerate 

Case-ll  (a) 


GDL-CL 

Interface 


Membrane-CL 

Interface 


Reference 

Agglomerate 


Case-ll  (b) 


Agglomerate 

Case-Ill  (a) 


GDL-CL 

Interface 


Membrane-CL 

Interface 


Reference 

Agglomerate 


Case-Ill  (b) 


Fig.  5.  Schematic  of  the  agglomerate  arrangements  in  the  cathode  catalyst  layer  in  part  (a)  and  the  computational  domain  in  part  (b).  Case-I  represents  in-line 
agglomerate  arrangement,  Case-II  represents  staggered  arrangement  in  y-direction,  and  Case-Ill  depicts  staggered  arrangement  in  both  y-  and  z-directions. 


is  observed  that  all  three  grid  levels  show  identical  results. 
A  close  inspection  of  these  results  show  about  0.5%  varia¬ 
tion  in  the  oxygen  concentration  and  about  0.2%  variation  in 
the  activation  overpotential  between  the  highest  and  the  low¬ 
est  level  meshes.  Nevertheless,  results  shown  in  Fig.  6  are 
almost  independent  of  the  grid  sizes.  Since  the  Case-II  and 
Case-Ill  consist  more  spherical  surfaces  than  the  Case-I  and 
also  for  better  accuracy,  results  are  presented  for  the  highest 
grid  level  (approximately  1  million  DOF)  in  all  the  subsequent 
figures. 


4.  Results  and  discussion 

The  accuracy  of  the  three-dimensional  finite  element  cal¬ 
culations  has  already  been  shown  in  the  previous  section 
by  comparing  three-dimensional  model  results  with  a  two- 
dimensional  axi-symmetric  model,  as  well  as  results  were 
compared  for  the  grid  sensitivity  for  the  three  different  grid 
levels.  In  this  section,  results  of  a  parametric  study  with  three- 
different  catalyst  agglomerate  arrangements  are  presented.  As 
shown  in  Fig.  5,  in  Case-I,  the  catalyst  agglomerates  are  con- 
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Fig.  6.  Variation  of  oxygen  concentration  in  the  cathode  catalyst  layer  along  x-axis  for  Case-I  in  a  PEM  fuel  cell  operating  at  80  °C  and  3  atm  in  part  (a)  and  the 
corresponding  variation  of  the  activation  overpotential  in  part  (b)  for  two  current  density  values  (J&  =  0.4  A  cm-2  and  0.6  A  cm-2).  Lines  represent  the  results  at  the 
lowest  degrees  of  freedom  (DOF),  whereas  the  symbols  represent  the  results  of  two  highest  degrees  of  freedom  as  indicated  in  the  figure. 


sidered  as  in-line  arrangement  in  all  directions.  For  Case-II 
and  Case-Ill,  two  staggered  arrangements  are  considered  as 
shown  in  Fig.  5.  Here,  staggered  arrangement  is  considered  in 
y-direction  for  Case-II,  and  for  Case-Ill,  both  y-  and  z-directions 
have  staggered  arrangements  of  the  catalyst  agglomerates.  For 
all  the  cases,  x-direction  is  considered  as  in-line  arrangement 
to  keep  symmetry  between  these  cases.  Due  to  the  symmetry 
in  x-direction,  all  the  results  are  presented  along  the  centerline 
between  the  two  reference  agglomerates  (lied  on  the  x-axis)  as 
shown  in  part  (b)  of  Fig.  5. 

4.1.  Case-I  ( in-line  arrangement ) 

Fig.  7  shows  the  profiles  of  oxygen  concentration  with  dif¬ 
ferent  current  density  values  for  Case-I.  Each  of  these  profiles  is 
plotted  inside  the  reference  agglomerates  along  the  x-axis  that 
is  identical  of  line  00  as  shown  in  Fig.  2.  Five  different  cur¬ 
rent  density  values  were  considered  as  indicated  in  the  legend. 


In  both  figures,  the  symbols  represent  the  oxygen  concentration 
profile  along  a  line  parallel  to  the  x-axis  at  y  =  Oandz  =  2.5  p,m 
(equal  to  agglomerate  radius)  for  J$  =  0. 1  A  cm-2,  i.e.,  along 
the  interface  between  two  agglomerates  on  xz-plane.  The  sim¬ 
ulation  parameters  used  to  estimate  the  oxygen  concentration 
are  listed  in  Table  1 .  Two  parts  of  this  figure  depict  two  differ¬ 
ent  combinations  of  operating  parameters,  namely,  T  =  50  °C 
and  P  =  1  atm,  and  T  —  80  °C  and  P  =  3  atm.  Here  x  =  0  rep¬ 
resents  the  membrane-CL  interface  and  x  =  10  fxm  represents 
the  GDL-CL  interface.  It  is  observed  that  the  oxygen  concentra¬ 
tion  decreases  along  the  centerline  from  the  GDL-CL  interface 
to  the  membrane-CL  interface.  For  low  current  densities,  the 
variation  in  the  concentration  is  less  whereas  for  high  current 
densities  an  oscillatory  behavior  is  observed  in  the  profile.  As 
expected,  the  minimum  oxygen  concentration  is  observed  at  the 
center  point  of  the  agglomerate  and  the  two  undulations  in  the 
profiles  represent  the  two  agglomerates.  Although  the  concen¬ 
tration  profile  shows  a  decreasing  behavior,  significant  amount 


Fig.  7.  Oxygen  concentration  profile  in  the  cathode  catalyst  layer  along  the  center  line  (x-axis  in  Fig.  5)  of  the  agglomerates  for  Case-I  in  a  PEM  fuel  cell  operating 
at  (a)  T  =  50  °C  and  P  =  1  atm,  and  (b)  T  =  80  °C  and  P  =  3  atm.  Each  line  represents  result  of  different  current  density  values  as  indicated  in  the  legend,  while 
the  symbols  show  the  oxygen  profile  along  a  line  parallel  to  x-axis  at  y  =  0  and  z  =  2.5  |xm  for  J§  =  0.1  A  cm-2.  All  other  parameters  are  listed  in  Table  1. 
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Fig.  8.  Distribution  of  the  activation  overpotential  in  the  cathode  catalyst  layer  for  Case-I  corresponding  to  the  oxygen  concentration  shown  in  Fig.  7.  Each  line 
represents  different  current  density  values  as  indicated  in  the  legend  while  part  (a)  for  T  =  50  °C  and  P  =  1  atm,  and  part  (b)  for  T  =  80  °C  and  P  =  3  atm. 


of  oxygen  concentration  still  exists  in  the  membrane-CL  inter¬ 
face  due  to  the  fast  oxygen  diffusion  through  the  void  region 
around  the  catalyst  agglomerates.  Further,  in  this  study,  an  ideal 
case  scenario  is  considered  when  there  is  no  flooding  outside 
the  catalyst  agglomerates,  whereas  the  catalyst  agglomerates 
are  partially  flooded.  Hence,  oxygen  diffusion  across  the  dry 
void  region  dominates  over  the  diffusion  through  the  partially 
flooded  catalyst  agglomerates  as  shown  by  symbols  in  both  fig¬ 
ure  for  J$  =  0.1  A  cm-2.  Here,  oxygen  concentration  is  almost 
constant  in  the  void  region  along  x-axis  due  to  the  favorable  oxy¬ 
gen  transport.  Only  variation  is  observed  at  the  contact  surfaces 
between  the  reference  and  surrounding  agglomerates.  Since  no 
variation  is  observed  in  the  oxygen  concentration  profile,  in  sub¬ 
sequent  figures,  results  in  the  void  region  have  not  been  reported. 
It  is  also  clear  from  Fig.  7  that  in  the  agglomerate,  oxygen  is 
transported  in  two  ways:  first  oxygen  diffuses  along  the  axial 
direction  or  the  thickness  of  the  catalyst  layer  in  the  void  region, 
then  from  the  void  region,  oxygen  diffuses  in  the  radial  direction 
towards  the  center  of  the  each  agglomerate. 

In  Fig.  8,  the  cathode  activation  overpotential  is  plotted 
as  a  function  of  spatial  coordinate  x  in  the  catalyst  layer  for 
different  current  density  values  as  indicated  in  the  legend.  Sim¬ 
ilar  to  Fig.  7,  two  parts  show  two  different  combinations  of 
the  operating  parameters  as  indicated  in  the  figure.  Each  line 
corresponds  to  the  activation  overpotential  for  the  oxygen  con¬ 
centration  shown  in  Fig.  7.  Identical  to  the  concentration  profile, 
the  variation  of  the  activation  overpotential  in  the  catalyst  layer 
is  small  at  low  current  densities,  whereas  the  activation  overpo¬ 
tential  decreases  rapidly  from  the  membrane-CL  interface  to  the 
GDL-CL  interface  for  the  higher  current  densities.  For  all  the 
current  densities,  the  activation  overpotential  is  higher  for  lower 
pressure  and  temperature  than  the  overpotential  for  the  higher 
pressure  and  temperature.  This  is  mainly  due  to  the  fast  oxygen 
transport  in  the  catalyst  layer  at  80  °C  and  3  atm. 

Fig.  9  illustrates  the  variations  of  the  reaction  rate  in  the  cat¬ 
alyst  layer  corresponding  to  the  oxygen  concentration  and  the 
activation  overpotential  shown  in  Figs.  7  and  8,  respectively, 
for  five  different  current  densities.  Here,  the  lines  represent 
the  results  for  T  =  50  °C  and  P  =  1  atm,  and  the  symbols  are 


T  =  80  °C  and  P  =  3  atm.  Surprisingly,  changing  the  operating 
condition  does  not  show  any  significant  effect  on  the  reac¬ 
tion  rate.  However,  slightly  higher  reaction  rate  is  observed 
for  T  =  80  °C  and  P  =  3  atm  at  the  center  of  the  agglomer¬ 
ates  for  high  current  densities.  These  similarities  show  that  the 
rate  of  the  electro-chemical  reaction  is  not  responsible  for  the 
difference  observed  in  the  activation  overpotential  in  Fig.  8  for 
different  operating  conditions,  which  is  solely  due  to  the  vari¬ 
ation  of  oxygen  concentration  in  the  catalyst  layer.  Although 
the  reactions  are  faster  at  higher  temperatures  and  pressures, 
here  it  has  not  been  significantly  observed  since  a  higher  operat¬ 
ing  temperature  and  pressure  is  known  to  reduce  the  activation 
overpotential  which  is  the  driving  force  for  the  electro-chemical 
reactions  occurring  in  the  fuel  cells. 

4. 2.  Case-II  ( staggered  arrangement-I) 

For  the  staggered  arrangements  of  catalyst  agglomerates, 
two  cases  were  considered.  In  the  first  staggered  arrangement, 
agglomerates  were  considered  as  staggered  in  y-direction,  i.e., 


x[gm] 

Fig.  9.  Variation  of  the  reaction  rate  in  the  cathode  catalyst  layer  along  x-axis 
for  Case-I.  Lines  represent  the  results  for  operating  conditions  of  T  =  50  °C  and 
P  =  1  atm,  and  the  symbols  depict  the  corresponding  results  for  T  =  80  °C  and 
P  =  3  atm. 
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Fig.  10.  Oxygen  concentration  profile  in  the  cathode  catalyst  layer  along  x-axis  for  Case-II  in  a  PEM  fuel  cell  operating  at:  (a)  T  =  50  °C  and  P  =  1  atm,  and  (b) 
T  =  80  °C  and  P  =  3  atm.  Each  line  represents  result  of  different  current  density  values  as  indicated  in  the  legend.  All  other  parameters  are  same  as  Case-I. 


uni-directional  staggered  arrangement,  as  shown  in  Fig.  5  as 
Case-II.  To  maintain  a  similarity  with  Case-I,  the  thickness  of 
the  catalyst  layer  is  kept  10  p,m.  Here,  the  reference  agglom¬ 
erates  (lied  on  the  x-axis)  are  considered  as  spherical,  whereas 
the  surrounding  agglomerates  can  be  either  spherical  or  hemi¬ 
spherical  to  maintain  the  thickness  of  the  catalyst  layer  same 
for  all  cases.  All  the  model  results  for  this  case  (Case-II)  are 
also  presented  along  the  centerline  of  the  two  middle  agglom¬ 
erates,  i.e.,  along  the  x-axis.  The  oxygen  concentration  profile 
in  the  catalyst  layer  for  Case-II  is  shown  with  different  cur¬ 
rent  density  values  in  Fig.  10.  All  the  simulation  parameters 
are  identical  of  Case-I.  Results  of  two  different  combination 
of  operating  parameters,  namely,  T  =  50  °C  and  P  =  1  atm, 
and  T  =  80  °C  and  P  =  3  atm  are  shown  in  Fig.  10  a  and  b, 
respectively.  Although  the  oxygen  concentration  at  the  inter¬ 
face  of  GDL-CL  for  similar  temperature  and  pressure  are  equal 
for  both  Case-I  and  Case-II;  for  Case-II,  a  smaller  value  of  oxy¬ 
gen  concentration  is  observed  at  the  interface  of  membrane-CL. 


This  is  reasonable,  since  in  Case-II,  agglomerates  are  staggered 
in  y-direction.  Hence,  Case-II  has  less  void  space  around  the 
agglomerates  compared  to  Case-I,  which  eventually  prevents 
faster  oxygen  diffusion  through  the  constricted  void  spaces.  Fur¬ 
thermore,  the  undulatory  profile  in  the  oxygen  concentration  is 
more  prominent  in  this  case.  Qualitatively,  oxygen  concentration 
profile  for  both  the  pressure  and  temperature  combinations  show 
similar  behavior  except  their  magnitudes.  Further  inspection  on 
the  values  of  oxygen  concentration  at  the  GDL-CL  interface 
shows  that  at  higher  temperature  and  pressure,  concentration  is 
higher  than  the  smaller  temperature  and  pressure  combination. 
This  is  mainly  due  to  the  faster  transport  processes  through  the 
un-flooded  GDL  at  higher  pressure  and  temperature,  though  oxy¬ 
gen  concentration  in  the  gas  channel  is  less  for  T  =  80  °C  and 
P  =  3  atm  due  to  the  higher  fraction  of  water  vapor.  Irrespec¬ 
tive  to  the  magnitude  of  the  oxygen  concentration  in  the  catalyst 
layer,  temperature  and  pressure  does  not  show  any  significant 
effect  on  the  profile  of  oxygen  concentration  in  the  catalyst  layer. 
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Fig.  11.  Distribution  of  the  activation  overpotential  in  the  cathode  catalyst  layer  for  Case-II  corresponding  to  the  oxygen  concentration  shown  in  Fig.  10.  Each  line 
represents  different  current  density  values  as  indicated  in  the  legend  while  part  (a)  for  T  =  50  °C  and  P  =  1  atm,  and  part  (b)  for  T  =  80  °C  and  P  =  3  atm. 
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Fig.  12.  Variation  of  the  reaction  rate  in  the  cathode  catalyst  layer  along  x-axis 
for  Case-II.  Lines  represent  the  results  for  operating  conditions  of  T  =  50°C  and 
P  =  1  atm,  and  the  symbols  represent  the  corresponding  results  for  T  =  80  °C 
and  P  =  3  atm. 

The  variations  observed  here  is  mainly  due  to  the  catalyst  layer 
structures,  or  in  the  other  words,  operating  conditions  dictate  the 
quantity  in  the  transport  process  whereas  agglomerate  structures 
dictate  the  quality  of  the  diffusion  in  the  transport  processes. 

Similar  to  the  Case-I,  the  results  of  Case-II  for  T  =  50  °C 
and  P  —  1  atm  show  higher  activation  overpotential  than  the 
corresponding  overpotential  for  T  =  80  °C  and  P  =  3  atm  as 
shown  in  Fig.  11.  Comparing  Fig.  1 1  with  Fig.  8  reveals  higher 
activation  loss  for  staggered  agglomerate  arrangements,  since 
lesser  path  available  for  the  oxygen  transport  due  to  the  stag¬ 
gered  agglomerate  orientation  in  the  catalyst  layer.  The  variation 
of  the  reaction  rate  in  the  cathode  catalyst  layer  for  the  Case- 
II  is  shown  in  Fig.  12.  Here,  the  lines  represent  the  results  for 
T  =  50  °C  and  P  —  1  atm,  and  the  symbols  represent  results 
corresponding  to  T  =  80  °C  and  P  =  3  atm  for  five  current 
densities  as  indicated  in  the  legend.  Here  a  distinct  difference  is 
observed  in  the  reaction  rate  profile  in  two  operating  conditions. 
For  higher  temperature  and  pressure,  the  rate  of  reaction  is  higher 
nearby  the  center  of  the  agglomerates  at  high  current  density  val¬ 
ues.  An  undulatory  nature  has  been  observed  in  the  reaction  rate 
profile,  while  a  crest  exists  at  the  interface  between  the  two  refer¬ 
ence  agglomerates.  This  undulatory  behavior  is  more  prominent 
in  higher  operating  parameter  values  and  high  current  densi¬ 
ties.  Furthermore,  in  some  instance,  the  higher  current  density 
shows  lower  reaction  rate  at  the  center  of  the  agglomerates  for 
T  =  50  °C  and  P  =  1  atm,  which  might  be  due  to  insufficient 
oxygen  available  on  the  surface  of  the  agglomerates  or  slow  dif¬ 
fusion  towards  the  center  of  the  agglomerates  at  higher  current 
densities. 

4. 3.  Case -III  ( staggered  arrangement-II) 

In  the  second  type  of  staggered  arrangement,  the  arrange¬ 
ments  of  catalyst  agglomerates  were  considered  as  staggered 


in  y-  and  z-directions,  i.e.,  bi-directional  staggered  arrange¬ 
ment,  as  shown  in  Fig.  5  as  Case-Ill.  Like  the  other  two  cases, 
all  the  results  are  presented  along  the  centerline  of  the  two 
reference  agglomerates  that  is  along  x-  axis  (see  Case-Ill  in 
Fig.  5).  Fig.  13  shows  the  variations  of  the  oxygen  concentration 
with  different  current  density  values  as  indicated  in  the  legend 
for  T  =  50  °C  and  P  =  1  atm  (Fig.  13a),  and  T  =  80  °C  and 
P  =  3  atm  (Fig.  13b).  The  oxygen  concentration  profiles  show 
almost  identical  behavior  like  Case-II.  However,  at  the  center 
of  the  agglomerates,  oxygen  concentration  is  slightly  higher 
than  the  Case-II.  Since  the  Case-Ill  has  staggered  structures 
of  agglomerate  arrangements  in  two  directions,  the  diffusion 
around  the  agglomerates  is  non-uniform.  This  non-uniformity 
eventually  provides  favorable  environment  for  the  oxygen  dif¬ 
fusion  in  the  radial  direction  of  the  agglomerates  for  Case-Ill. 
Furthermore,  better  oxygen  diffusion  also  reduces  the  activa¬ 
tion  overpotential  for  the  Case-Ill  as  shown  in  Fig.  14  for  both 
operating  conditions  compared  to  the  Case-II.  This  can  be  bet¬ 
ter  explained  by  visualizing  the  agglomerate  arrangements.  For 
Case-II,  staggered  arrangement  exists  in  one  direction,  there¬ 
fore,  oxygen  diffusion  outside  the  agglomerates  is  faster  in  the 
direction  where  in-line  arrangements  exist,  and  slower  in  the 
staggered  direction.  When  it  comes  to  the  agglomerates  sur¬ 
face,  then  more  oxygen  is  available  in  certain  areas,  whereas 
not  all  oxygen  can  diffuse  inside  the  agglomerate  through  the 
radial  direction.  This  can  be  visualized  by  plotting  oxygen  con¬ 
centration  in  the  void  region  for  Case-II  and  Case-Ill.  Fig.  15 
shows  oxygen  concentration  profile  outside  the  agglomerates 
along  a  line  parallel  to  the  x-axis  at  y  =  0  and  z  =  2.5  |mm  for 
T  =  80  °C  and  P  =  3  atm.  Here,  the  lines  depict  results  for 
Case-II  and  the  symbols  for  Case-Ill  for  different  current  den¬ 
sity  values  as  indicated  in  the  legend.  As  observed  in  Fig.  15, 
oxygen  concentration  on  the  catalyst  surface  for  Case-Ill  is 
higher  than  Case-II  at  the  interfaces  between  two  agglomerates 
(2  [xm  <  x  <3  [Jim  and  7  p,m  <  x  <8  p,m)  for  all  the  current 
density  values.  Therefore,  more  oxygen  available  for  Case-Ill 
to  diffuse  in  the  y-direction  that  is  due  to  the  staggered  arrange¬ 
ments  in  the  z-direction.  Furthermore,  due  to  the  higher  oxygen 
concentration  in  certain  areas,  diffusion  will  be  faster  for  Case- 
III  in  y-direction.  Also,  due  to  the  in-line  arrangements  in  two 
directions,  less  oxygen  available  in  certain  areas  on  the  agglom¬ 
erate  surface  for  Case-II  whereas  capacity  of  diffusion  through 
the  radial  direction  is  more.  When  it  comes  to  the  diffusion 
through  the  individual  agglomerate,  all  the  geometries  have 
same  composition  inside  the  agglomerate;  hence,  the  entire  dif¬ 
fusion  processes  is  control  by  the  amount  of  oxygen  available 
at  the  surface  of  the  agglomerates  and  how  the  concentration 
is  distributed  over  the  agglomerate  surfaces.  In  other  words, 
capacity  of  the  diffusion  is  higher  than  the  amount  of  oxygen 
available  on  the  agglomerate  surfaces  for  Case-II.  Combining 
these  effects  eventually  lower  the  oxygen  concentration  at  the 
agglomerate  center  in  Case-II  (see  Fig.  10).  For  Case-Ill,  since 
staggered  structures  exist  in  two  directions,  available  oxygen  on 
the  surface  of  the  agglomerates  is  higher  in  the  y-direction  than 
Case-II.  This  eventually  enhances  the  diffusion  process  in  the 
radial  direction  due  to  the  higher  oxygen  availability  and  low¬ 
ers  the  activation  overpotential  as  shown  in  Fig.  14.  It  might 
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Fig.  13.  Oxygen  concentration  profile  in  the  cathode  catalyst  layer  along  the  x-axis  for  Case-Ill  in  a  PEM  fuel  cell  operating  at:  (a)  T  =  50  °C  and  P  =  1  atm,  and 
(b)  T  =  80  °C  and  P  =  3  atm.  Each  line  represents  result  of  different  current  density  values  as  indicated  in  the  legend. 
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Fig.  14.  Distribution  of  the  activation  overpotential  in  the  cathode  catalyst  layer  for  Case-Ill  corresponding  to  the  oxygen  concentration  shown  in  Fig.  13.  Each  line 
represents  different  current  density  values  as  indicated  in  the  legend  while  part  (a)  for  T  =  50  °C  and  P  =  1  atm,  and  part  (b)  for  T  =  80  °C  and  P  =  3  atm. 


be  questionable,  though  Case-Ill  has  higher  non-uniformity  in 
the  catalyst  arrangements  then  why  Case-Ill  show  better  dif¬ 
fusion  results.  The  simple  answer  will  be  it  provides  better 
passage  in  the  void  region  in  z-direction  due  to  the  staggered 
arrangements. 

In  the  reaction  rate  of  Case-Ill,  similar  behavior  is  observed 
like  Case-II.  The  plot  of  the  reaction  rate  for  Case-Ill  is  shown  in 
Fig.  16  for  T  =  50°CandF  =  1  atm  as  lines,  and  for  T  =  80  °C 
and  P  =  3  atm  as  symbols.  For  all  the  current  densities,  the  reac¬ 
tion  rate  is  higher  for  higher  operating  parameters  and  lower  for 
lower  operating  parameters.  The  local  variation  in  the  reaction 
rate  observed  here  mainly  due  to  the  local  change  of  oxygen 
concentration  observed  in  Fig.  13  for  different  operating  con¬ 
ditions.  Like  the  activation  overpotential,  the  rate  of  reaction  is 
also  lower  for  Case-Ill  compared  to  Case-II  as  shown  in  Fig.  17. 
In  Fig.  17,  the  reaction  rates  are  plotted  for  all  the  three  cases 
along  the  v-axis  for  J§  =  0.6  A  cm-2  for  the  fuel  cell  oper¬ 
ating  at  T  =  80  °C  and  P  =  3  atm.  For  the  entire  thickness  of 
the  catalyst  layer,  reaction  rate  is  highest  for  the  Case-II  and 


x[pm] 

Fig.  15.  Oxygen  concentration  profile  in  the  cathode  catalyst  layer  along  a  line 
parallel  to  x-axis  at  y  =  0  and  z  =  2.5  [xm  for  T  =  80  °C  and  P  =  3  atm.  Lines 
represent  the  result  for  Case-II,  while  the  symbols  represent  Case-Ill  for  five 
different  current  density  values  as  indicated  in  the  legend. 
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Fig.  16.  Variation  of  the  reaction  rate  in  the  cathode  catalyst  layer  along  x-axis  for 
Case-Ill.  Lines  represent  the  results  for  operating  conditions  of  T  =  50  °C  and 
P  =  1  atm,  and  the  symbols  represent  the  results  for  T  =  80  °C  and  P  =  3  atm. 

lowest  for  the  Case-I.  This  might  be  another  possible  cause  for 
the  lower  oxygen  concentration  observed  at  the  center  of  the 
agglomerates  for  Case-II  (cf.  Fig.  10)  compared  to  Case-Ill. 
In  addition,  concentration,  reaction  rate,  and  activation  over¬ 
potential  are  coupled;  therefore,  it  is  required  to  optimize  those 
variables  in  order  to  find  better  cell  performance  and  design. 
Nevertheless,  the  results  obtained  from  the  present  investigation 
reveals  a  considerable  insight  on  how  the  governing  parameters 
for  PEM  fuel  cells  change  with  the  structures  of  the  catalyst 
layer  as  well  as  with  the  operating  conditions.  Like  the  Case-II 
shows  highest  reaction  rate  that  means  the  speed  of  the  chemical 
reaction  is  faster,  hence  less  catalyst  will  be  required  to  pro¬ 
mote  the  electro-chemical  reaction.  On  the  other  hand,  Case-II 
also  shows  higher  activation  losses.  Therefore,  the  results  pre¬ 
sented  in  this  paper  providing  information,  and  the  direction 
when  and  where  the  optimization  is  possible  or  at  what  extent  it  is 
possible. 


x[nm] 

Fig.  17.  Comparison  between  the  reaction  rates  at  J$  =  0.6  Acm_2for  the  fuel 
cell  operating  at  T  =  80  °C  and  P  =  3  atm.  Each  line  represents  the  reaction 
rate  along  the  x-axis  for  three  cases  as  indicated  in  the  legend. 


5.  Conclusions 

In  this  paper,  a  three-dimensional  agglomerate  model  for 
the  cathode  catalyst  layer  in  a  PEM  fuel  cell  has  been  devel¬ 
oped  to  study  the  effect  of  catalyst  agglomerate  arrangements. 
A  comprehensive  validation  has  also  been  shown  that  indi¬ 
cates  complex  staggered  arrangement  of  catalyst  agglomerates 
in  three-dimensional  catalyst  layer  modeling  is  possible.  Consid¬ 
ering  the  three  different  arrangements  of  catalyst  agglomerates 
investigated  here,  it  is  found  that  the  agglomerate  arrangements 
in  the  catalyst  layer  have  significant  effect  on  the  oxygen  trans¬ 
port  processes  in  the  catalyst  layer,  consequently  the  activation 
losses  and  the  rate  of  reactions.  Therefore,  in  the  fuel  cell  mod¬ 
eling,  detailed  catalyst  layer  structures  should  be  considered 
to  capture  the  true  nature  of  the  PEM  fuel  cell  performance, 
particularly,  its  activation  losses.  It  has  also  been  observed 
that  in-line  arrangement  of  the  catalyst  agglomerates  yields 
lowest  activation  loss,  although  in  reality  catalyst  agglomer¬ 
ates  may  be  staggered  in  all  directions.  For  example,  almost 
17%  of  reduction  in  the  activation  overpotential  is  possible  by 
using  in-line  arrangement  of  the  catalyst  agglomerates  at  the 
current  density  of  0.4  A  cm-2  for  T  =  80  °C  and  P  =  3  atm. 
Among  the  two  staggered  arrangements,  bi-directional  (Case- 
III)  staggered  arrangement  shows  less  activation  overpotential 
than  uni-directional  (Case-II)  staggered  arrangement.  In  fact, 
3-4%  of  activation  losses  can  be  recovered  (considering  Case-II 
with  Case-Ill)  just  by  changing  the  agglomerate  arrangements. 
Hence,  it  might  be  possible  to  identify  an  optimum  catalyst  layer 
structures  by  further  investigating  such  staggered  arrangements. 
Furthermore,  it  seems  that  the  present  catalyst  layer  models  will 
be  useful  in  optimizing  the  catalyst  layer  structures  and  catalyst 
distribution. 
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